Identification of key ferroptosis genes in diabetic retinopathy based on bioinformatics analysis

Objectives Diabetic retinopathy (DR) is a retinal microvascular disease associated with diabetes. Ferroptosis is a new type of programmed cell death that may participate in the occurrence and development of DR. Therefore, this study aimed to identify the DR ferroptosis-related genes by bioinformatics methods. Methods The RNAseq data of DR and healthy control retinas were downloaded from the gene expression synthesis (GEO) database and analyzed using the R package DESeq2. The key modules were obtained using the WGCNA algorithm, and their genes were intersected with ferroptosis-related genes in the FerrDb database to obtain differentially expressed ferroptosis-related genes (DE-FRGs). Enrichment analysis was conducted to understand the function and enrichment pathways of ferroptosis genes in DR, and hub genes were identified by protein-protein interaction (PPI) analysis. The diagnostic accuracy of hub genes for DR was evaluated according to the area under the ROC curve. The TRRUST database was then used to predict the regulatory relationship between transcription factors and target genes, with the mirDIP, ENCORI, RNAnter, RNA22, miRWalk and miRDB databases used to predict the regulatory relationship between miRNAs and target genes. Finally, another data set was used to verify the hub genes. Results In total, 52 ferroptosis-related DEGs (43 up-regulated and 9 down-regulated) were identified using 15 DR samples and 3 control samples and were shown to be significantly enriched in the intrinsic apoptotic signaling pathway, autophagosome, iron ion binding and p53 signaling pathway. Seven hub genes of DR ferroptosis were identified through PPI network analysis, but only HMOX1 and PTGS2 were differentially expressed in another data set. The miRNAs prediction showed that hsa-miR-873-5p was the key miRNA regulating HMOX1, while hsa-miR-624-5p and hsa-miR-542-3p were the key miRNAs regulating PTGS2. Furthermore, HMOX1 and PTGS2 were regulated by 13 and 20 transcription factors, respectively. Conclusion The hub genes HMOX1 and PTGS2, and their associated transcription factors and miRNAs, may be involved in ferroptosis in diabetic retinopathy. Therefore, the specific mechanism is worthy of further investigation.


Introduction
Diabetic retinopathy (DR) is a retinal microvascular disease associated with diabetes, accounting for an estimated 4.8% of global blindness, and has become the main cause of acquired visual loss in middle-aged people worldwide [1]. DR is classified as non-proliferative or proliferative DR according to the development stage [2]. Oxidative stress and inflammation play an important role in DR pathophysiology [3] but the pathogenesis and mechanism of DR remain unclear.
Ferroptosis is a new mode of programmed cell death that is different from apoptosis and involves high iron-dependent lipid peroxidation [4]. The cumulation of lipid peroxides leads to the loss of selective cell membrane permeability [5]. However, glutathione peroxidase 4 (GPX4) plays an important role in protecting cell membrane from peroxidation damage that can prevent the cell membrane from being affected by oxidation [6]. When ferroptosis occurs, the antioxidant glutathione is exhausted, leading to GPX4 failure and ultimately, fatal accumulation of lipid peroxides [7].
The high glucose environment can inhibit the growth of human retinal capillary endothelial cells, and ferroptosis can increase this effect, possibly related to GPX4 ubiquitination promoted by TRIM46 [8]. The fatty acid binding protein 4 (FABP4) can inhibit lipid peroxidation and oxidative stress by regulating ferroptosis, thereby reducing retinal damage in DR [9]. In addition, non-coding RNA (ncRNA), such as circular RNA and miRNA, also participate in the ferroptosis of DR [10,11]; thus, ferroptosis is closely related to DR initiation and progression, but the underlying mechanism of ferroptosis in DR remains unknown.
Therefore, this study identified differentially expressed genes (DEGs) in the retina of DR patients and normal retinas and then matched these DEGs to the ferroptosis dataset to acquire differentially expressed ferroptosis-related genes (DE-FRGs) to determine the mechanism of DR development. Our findings reveal the role of ferroptosis in DR, which may become a target of clinical drug therapy in the future.

Data collection and acquisition of ferroptosis-related gene
GEO (https://www.ncbi.nlm.nih.gov/geo/) belongs to public databases. Users can download relevant data for free for research and publish relevant articles [12]. In this study, the GSE102485 mRNA expression profile dataset was downloaded from GEO. GSE102485 contains 19 specimens of diabetic retinopathy and 3 normal retinas and is based on the GPL18573 platform (Illumina NextSeq 500, Homo sapiens). The public FerrDb (http://www.zhounan. org/ferrdb) database for ferroptosis-related genes (FRGs) that promote, inhibit, or mark ferroptosis was also searched, and after the removal of duplicate genes, 232 FRGs were finally obtained for subsequent analysis.

Identification of DEGs
R software (version 4.1.3; https://www.r-project.org/) and the Bioconductor software package (http://www.bioconductor.org/)) were used to correct and analyze the original data. RNA-seq data processing and normalization were performed using DESeq2 packages. Principal component analysis (PCA) was used to verify the repeatability of the GSE102485 data. The standard of statistical significance was |log2FC| > 1 and adjusted p < 0.05. DEGs were visualized in a volcano map based on the "ggplot2" software package.

Weighted gene co-expression network analysis
Based on the scale-free topology criterion, a co-expression network of DEGs was constructed. First, all thresholds were analyzed using the WGCNA software package in R software to determine the soft threshold power. Then, a weighted co-expression network was constructed, and the classifier was clustered into multiple modules with different colored labels to determine the correlation between each module and the control. The module most related to diabetes was considered the key module for further analysis.

Differential expression analysis of ferroptosis-related genes
The differentially expressed ferroptosis-related genes (DE-FRGs) were cross genes between key module genes and FRGs. At the same time, determine whether DE-FRGs are genes of single gene symbols. We used the "Venn Diagram" package of R software to draw Venn diagram to show the number of DE-FRGs. Used the "ggplot2" package to display the expression of DE-FRGs in the heat-map.

GO and KEGG analysis of DE-FGRs
Functional annotations and pathway enrichment for GO biological processes and KEGG annotation were performed using the "Cluster analyzer" software package. The enrichment results were sorted according to the adjusted P value. The enrichment bar chart shown the first 10 results.

Construction of a PPI network of DE-FRGs and hub gene identification
The PPI network was constructed using the search tool of STRING online database (https:// cn.string-db.org/) and visualized by Cytoscape software. The top 10 genes of the PPI network were determined as the hub genes, which were calculated based on the maximal clique centrality (MCC), maximum neighborhood component (MNC), degree, edge percolated component (EPC) and Closeness algorithms y utilizing the cytohubba plug-in. The expression of hub genes was visualized in a boxplot, and the area under the curve (AUC) was measured by ROC curve analysis of hub genes. The black line diagram was drawn using the "circlize" package of R software to show the hub genes correlations.

TFs-Genes-miRNAs interaction networks
Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) is the largest freely available database of human transcription factors (TFs) and genes interactions [13]. We obtained the TFs regulating hub gene from TRRUST database. Then mirDIP, ENCORI, RNAnter, RNA22, miRWalk and miRDB online databases were used to construct the miRNAs-genes cooperative regulation network. The TFs-genes-miRNAs network was visualized by Cytoscape software.

Verification of hub genes
The GSE94019 database from the GEO database was used to verify the accuracy of the identified hub genes. The "limma" package of R software was used to calibrate the data set, and differences between DR and controls were compared using a Wilcoxon test ( � p-value < 0.05 and �� p-value < 0.01). The box graph was drawn using R software.

Identification of DEGs
PCA was used to detect whether biological replicates in the same treatment group were clustered together and whether samples in different treatment groups were separated from each other. The analysis results showed that the GSE102485 data set had well repeatability ( Fig 1A). It was worth noting that PCA also showed that four DR samples were outliers, so they were excluded from the analysis. Subsequently, the genes of 15 DR retinal tissues and 3 normal retinal tissues were used for differential gene analysis. Used the adjusted p value < 0.05 and | log2FC| > 1 as criterion. In all, 4876 DEGs were identified, including 2979 up-regulated genes and 1897 down-regulated genes. The differential expression of these 4876 DEGs in the DR and the control groups is shown in the volcanic map ( Fig 1B).

WGCNA analysis
The WGCNA package in R software was used to further process 4876 identified DEGs. When building the sample tree, there were no abnormal samples, so no samples were removed ( Fig  2A). The soft threshold power of 18 was used to establish a scale-free co-expression network (scale-free R2 > 0.8) (Fig 2B). The clusters were divided into seven modules, blue, brown, green, red, yellow, turquoise and green, with a minimum module size of � 30 ( Fig 2C). The correlation between each module and DR was determined (Fig 2D), showing that blue (-0.88, P < 0.0001) and turquoise (0.95, P < 0.0001) were the most negative and positive modules related to DR, respectively. These two modules were the first two modules significantly related to clinical characteristics. The blue (cor = 0.7, P = 2.6e-151) and turquoise (cor = 0.91, P < 1e-200) modules showed a significant positive correlation between MM and GS of the target genes ( Fig 2E and 2F). Therefore, blue module and turquoise module were considered as key modules.

Identification of DE-FRGs
The data of 232 genes from the FerrDb database were crossed with the genes of key modules to determine DE-FRGs to identify. A total of 43 positively-related genes and 9 negative related genes were found. These 52 DE-FRGs are all single gene symbols, and the heat-map and Venn diagram are shown in Fig 3. These DE-FRGs were further classified as ferroptosis driver, suppressor, or marker genes (Table 1).

GO and KEGG analysis of DE-FRGs
DE-FRGs were enriched and analyzed by R software to define their underlying physiological functions. The enrichment histogram shows the top 10 results. The most significant enrichment terms were the intrinsic apoptotic signaling pathway, cellular response to external stimulus, response to oxidative stress, response to nutrient levels(biological process) (Fig 4A), secondary lysosome, the autolysosome, autophagosome, membrane raft (cellular component) (Fig 4B), iron ion binding, antioxidant activity, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding (molecular function)( Fig 4C). As expected, the KEGG pathway analysis showed that these DE-FRGs were significantly enriched in ferroptosis, the p53 signaling pathway, and endocrine resistance (Fig 4D).

PPI network of DE-FRGs
Protein-Protein interaction networks were created to demonstrate the interactions between DE-FRGs ( Fig 5A). We get a PPI network with 42 nodes and 115 edges. Ten of the 52 genes are not related to other molecules and do not form a molecular network. There are 47 up-regulated genes and 5 down-regulated genes in PPI network. Table 2

Hub gene detection
We discuss the diagnostic ability of these seven genes in different patients, and draw the ROC curve. The results show that the AUC of HMOX1, PTGS2, EGFR, CAV1, TLR4, MAPK8 and CDKN2A were 0.911, 1.000, 0.911, 0.978, 1.000, 1.000 and 0.889 respectively when distinguishing DR patients from normal controls (Fig 5B). The results show that 7 genes as new biomarkers have high diagnostic accuracy. The chord diagram shows that there is a strong correlation between the seven genes ( Fig 5C).

Verification of hub genes
The "ggpubr" package of R software was used to view the data distribution, and the values of each group centered on the median are shown in Fig 7A. This shows that the data in GSE94019 are uniform and comparable, so the samples passed the quality test. The comparison of the DEGs in the DR and the control groups showed that only two genes (HMOX1 and PTGS2) were different. The expression of hub genes in each group is shown in Fig 7B.

Discussion
Diabetic retinopathy (DR) is a special microvascular complication of diabetes and an important cause of blindness in many countries [14]. In the early stage, that is, non-proliferative DR,  the blood vessels in the retina are weakened, resulting in fluid leakage, macular swelling, and eventually blurred vision, whereas, in late proliferative DR, abnormal neovascularization on the surface of the retina seriously impairs vision [8]. Ii is predicted that by 2040, the number of DR patients will reach 191million globally, bringing great pressure to individuals, families, and society [15]. Therefore, elucidating the pathogenesis of DR is important to reduce the risk of blindness in diabetic patients.  Ferroptosis is a new mode of programmed cell death involving the accumulation of lipid peroxides that damage the cell membrane and is considered different from apoptosis, necrosis, and autophagy [16]. Li et al. used bioinformatics methods to identify the key genes of ferroptosis-related to Calcic aortic valve disease and proved in vitro that these hub genes are differentially expressed in Calcic aortic valve disease and normal people [17]. Previous studies have shown that ferroptosis is associated with aging retinopathy [18], light-induced retinal degeneration [19], retinitis pigmentosa [20], and glaucoma damage [21]. Therefore, we applied bioinformatics methods to explore hub genes involved in DR ferroptosis.
The present study showed that the expression of ferroptosis genes was different between the diabetic retinopathy and normal control groups. PCA indicated that the samples of DR group and control group were obviously different. There were 52 DE-FRGs that were involved in apoptosis, iron binding, autophagy, and other processes. The PPI network identified seven hub genes related to ferroptosis in DR, including HMOX1, PTGS2, EGFR, CAV1, TLR4, MAPK8, and CDKN2A. However, after verification with another data set, only HMOX1 and PTGS2 were shown to be differentially expressed between the DR and control groups. In addition, when testing the diagnostic accuracy of hub genes, the AUC of PTGS2 was equal to 1, possibly due to the small number of samples included in the study. The AUC of a subset of genes obtained using machine learning in the study by Li et al also equals 1 [22], which may indicate that an AUC of 1 is a relatively reliable diagnostic model in small-sample studies. In the future, further verification should be conducted using a larger sample size.
Heme oxygenase-1 (HMOX1)is the rate-limiting enzyme involved in the heme oxygenase reaction pathway, degrading heme into carbon monoxide (CO), ferrous (Fe 2+ ) and biliverdin IXα [23]. The increased HMOX1 in a mouse model reduced the neuronal cell death induced by oxidative stress [24]; however, excessive activation of HMOX1 can cause glioma cell death [25]. In retinopathy, the protective or damaging effect of HMOX1 on the retina depends on the expression of HMOX1 [26]. Our study shows that the high HMOX1 expression may be related to DR, but the role of HMOX1 in DR is still unclear. Prostaglandin endoperoxide synthase (PTGS) is a rate-limiting enzyme in the arachidonic acid synthesis of prostaglandins (PGs). It is expressed exists widely in mouse, rat and human retinas, including prostaglandin endoperoxide synthase-1 and prostaglandin endoperoxide synthase-2 (PTGS1 and PTGS2) [27]. PTGS2, which can be induced by cytokines, mitogens and endotoxins, is the immediate early genetic product of inflammation [28] and may be involved in DR by regulating the antiangiogenic factor TSP-1 and its receptor CD36 on endothelial cells [29]. Our study suggests that PTGS2 may be involved in the process of ferroptosis causing DR. Based on the close relationship between the high glucose environment and inflammation, PTGS2 and inflammation, the specific mechanism of PTGS2 participating in DR is worthy of further study.
MicroRNA (miRNA) is a small non-coding RNA involved in post-transcriptional gene regulation by degrading or inhibiting the translation of target genes [30]. TFs as modular proteins can bind to the DNA-binding domain in the promoter region of target genes to regulate transcription [31]. In our study, hsa-miR-873-5p was the key miRNA regulating HMOX1, while hsa-miR-624-5p and hsa-miR-542-3p were the key miRNAs regulating PTGS2. TFs prediction showed that HMOX1 and PTGS2 were regulated by 13 and 20 TFs respectively. It has been shown that miR-542-3p promotes the rapid degradation of PTGS2 mRNA, which partly supports our bioinformatics prediction [32]. MiR-624-5p [33][34][35] and miR-873-5p [36][37][38] have been studied in cancer, but their role in DR has not been investigated. According to the analysis of published literature, relevant TFs PPARG [39,40], STAT1 [41], AR [42], SP1 [43], STAT3 [44], RELA [45] and EGR1 [46] are closely related to DR. The above studies reflect that our predictions are relatively reliable. Our research suggests that TFs genes miRNAs connection may play a role in DR, but a lot of research is still needed to explore its specific mechanism.

Conclusion
In conclusion, we found that HMOX1 and PTGS2 are hub genes involved in ferroptosis process of diabetic retinopathy. Our study is helpful to increase the understanding of researchers on the pathogenesis of diabetic retinopathy. However, how these ferroptosis-related genes play a role needs to be explored by future in vivo and in vitro experiments.